home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / gauss2dfit.pro < prev    next >
Text File  |  1997-07-08  |  8KB  |  230 lines
  1. ; $Id: gauss2dfit.pro,v 1.3 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1995-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5. ;
  6.  
  7. PRO    GAUSS2_FUNCT, X, A, F, PDER
  8. ;+
  9. ; NAME:
  10. ;    GAUSS2_FUNCT
  11. ; PURPOSE:
  12. ;    Evaluate function for gauss2fit.
  13. ; CALLING SEQUENCE:
  14. ;    FUNCT,X,A,F,PDER
  15. ; INPUTS:
  16. ;    X = values of independent variables, encoded as: [nx, ny, x, y]
  17. ;    A = parameters of equation described below.
  18. ; OUTPUTS:
  19. ;    F = value of function at each X(i,j), Y(i,j).
  20. ;    Function is:
  21. ;        F(x,y) = A0 + A1*EXP(-U/2)
  22. ;        where: U= (yp/A2)^2 + (xp/A3)^2
  23. ;
  24. ;      If A has 7 elements a rotation of the ellipse is present and:
  25. ;        xp = (x-A4) * cos(A6) - (y-A5) * sin(A6)
  26. ;        yp = (x-A4) * sin(A6) + (y-A5) * cos(A6)
  27. ;      If A has 6 elements, A6 (theta) is 0, the major and minor axes
  28. ;      of the ellipse are parallel to the XY axes, and:
  29. ;        xp = (x-A4)   and   yp = (x-A5)
  30. ;
  31. ; Optional output parameters:
  32. ;    PDER = (n_elements(z),6 or 7) array containing the
  33. ;        partial derivatives.  pder(i,j) = derivative
  34. ;        at ith point w/respect to jth parameter.
  35. ; PROCEDURE:
  36. ;    Evaluate the function and then if requested, eval partials.
  37. ;
  38. ; MODIFICATION HISTORY:
  39. ;    WRITTEN, DMS, RSI, June, 1995.
  40. ;-
  41.  
  42. nx = long(x[0])        ;Retrieve X and Y vectors
  43. ny = long(x[1])
  44.  
  45. tilt = n_elements(a) eq 7    ;TRUE if angle present.
  46. if tilt then begin        ;Rotate?
  47.     xp = (x[2:nx+1]-a[4]) # replicate(1.0, ny)    ;Expand X values
  48.     yp = replicate(1.0, nx) # (x[nx+2:*]-a[5])    ;expand Y values
  49.     s = sin(A[6]) & c = cos(A[6])
  50.     t =  xp * (c/a[2]) - yp * (s/a[2])
  51.     yp = xp * (s/a[3]) + yp * (c/a[3])
  52.     xp = temporary(t)
  53. endif else begin
  54.     xp = (x[2:nx+1]-a[4]) # replicate(1.0/a[2], ny)    ;Expand X values
  55.     yp = replicate(1.0/a[3], nx) # (x[nx+2:*]-a[5])    ;expand Y values
  56.     s = 0.0 & c = 1.0
  57. endelse
  58.  
  59. n = nx * ny
  60. u = reform(exp(-0.5 * (xp^2 + yp^2)), n)    ;Exp() term, Make it 1D
  61. F = a[0] + a[1] * u
  62.  
  63. if n_params(0) le 3 then return ;need partial?  No.
  64.  
  65. PDER = FLTARR(n, n_elements(a))    ;YES, make partial array.
  66. PDER[*,0] = 1.0            ;And fill.
  67. pder[0,1] = u
  68. u = a[1] * u            ;Common term for the rest of the partials
  69. pder[0,2] = u * xp^2 / a[2]
  70. pder[0,3] = u * yp^2 / a[3]
  71. pder[0,4] = u * (c/a[2] * xp + s/a[3] * yp)
  72. pder[0,5] = u * (-s/a[2] * xp + c/a[3] * yp)
  73. if tilt then pder[0,6] = u * xp*yp*(a[2]/a[3]-a[3]/a[2])
  74. END
  75.  
  76.  
  77.  
  78. Function Gauss2dfit, z, a, x, y, NEGATIVE = neg, TILT=tilt
  79. ;+
  80. ; NAME:
  81. ;    GAUSS2DFIT
  82. ;
  83. ; PURPOSE:
  84. ;     Fit a 2 dimensional elliptical gaussian equation to rectilinearly
  85. ;    gridded data.
  86. ;        Z = F(x,y) where:
  87. ;         F(x,y) = A0 + A1*EXP(-U/2)
  88. ;       And the elliptical function is:
  89. ;        U= (x'/a)^2 + (y'/b)^2
  90. ;    The parameters of the ellipse U are:
  91. ;       Axis lengths are 2*a and 2*b, in the unrotated X and Y axes,
  92. ;        respectively.
  93. ;       Center is at (h,k).
  94. ;       Rotation of T radians from the X axis, in the CLOCKWISE direction.
  95. ;       The rotated coordinate system is defined as:
  96. ;        x' = (x-h) * cos(T) - (y-k) * sin(T)  <rotate by T about (h,k)>
  97. ;        y' = (x-h) * sin(T) + (y-k) * cos(T)
  98. ;
  99. ;    The rotation is optional, and may be forced to 0, making the major/
  100. ;    minor axes of the ellipse parallel to the X and Y axes.
  101. ;
  102. ;    The coefficients of the function, are returned in a seven
  103. ;    element vector:
  104. ;    a(0) = A0 = constant term.
  105. ;    a(1) = A1 = scale factor.
  106. ;    a(2) = a = width of gaussian in X direction.
  107. ;    a(3) = b = width of gaussian in Y direction.
  108. ;    a(4) = h = center X location.
  109. ;    a(5) = k = center Y location.
  110. ;    a(6) = T = Theta the rotation of the ellipse from the X axis
  111. ;        in radians, counterclockwise.
  112. ;
  113. ;
  114. ; CATEGORY:
  115. ;    curve / data fitting
  116. ;
  117. ; CALLING SEQUENCE:
  118. ;    Result = GAUSS2DFIT(z, a [,x,y])
  119. ;
  120. ; INPUTS:
  121. ;    Z = dependent variable in a 2D array dimensioned (Nx, Ny).  Gridding
  122. ;        must be rectilinear.
  123. ;    X = optional Nx element vector containing X values of Z.  X(i) = X value
  124. ;        for Z(i,j).  If omitted, a regular grid in X is assumed,
  125. ;        and the X location of Z(i,j) = i.
  126. ;    Y = optional Ny element vector containing Y values of Z.  Y(j) = Y value
  127. ;        for Z(i,j).  If omitted, a regular grid in Y is assumed,
  128. ;        and the Y location of Z(i,j) = j.
  129. ;
  130. ; Optional Keyword Parameters:
  131. ;    NEGATIVE = if set, implies that the gaussian to be fitted
  132. ;        is a valley (such as an absorption line).
  133. ;        By default, a peak is fit.
  134. ;    TILT = if set to  1, allow the orientation of the major/minor axes of 
  135. ;        the ellipse to be unrestricted.  The default is that
  136. ;        the axes of the ellipse must be parallel to the X-Y axes.
  137. ;        In this case, A(6) is always returned as 0.
  138. ;
  139. ; OUTPUTS:
  140. ;    The fitted function is returned.
  141. ; OUTPUT PARAMETERS:
  142. ;    A:    The coefficients of the fit.  A is a seven element vector as
  143. ;        described under PURPOSE.
  144. ;
  145. ; COMMON BLOCKS:
  146. ;    None.
  147. ; SIDE EFFECTS:
  148. ;    None.
  149. ; RESTRICTIONS:
  150. ;    Timing:  Approximately 4 seconds for a 128 x 128 array, on a 
  151. ;        Sun SPARC LX.  Time required is roughly proportional to the 
  152. ;        number of elements in Z.
  153. ;
  154. ; PROCEDURE:
  155. ;    The peak/valley is found by first smoothing Z and then finding the
  156. ;    maximum or minimum respectively.  Then GAUSSFIT is applied to the row
  157. ;    and column running through the peak/valley to estimate the parameters
  158. ;    of the Gaussian in X and Y.  Finally, CURVEFIT is used to fit the 2D
  159. ;    Gaussian to the data.
  160. ;
  161. ;    Be sure that the 2D array to be fit contains the entire Peak/Valley
  162. ;    out to at least 5 to 8 half-widths, or the curve-fitter may not
  163. ;    converge.
  164. ;
  165. ; EXAMPLE:  This example creates a 2D gaussian, adds random noise
  166. ;    and then applies GAUSS2DFIT:
  167. ;    nx = 128        ;Size of array
  168. ;    ny = 100
  169. ;    ;**  Offs Scale X width Y width X cen Y cen  **
  170. ;    ;**   A0  A1    a       b       h       k    **
  171. ;    a = [ 5., 10., nx/6.,  ny/10., nx/2., .6*ny]  ;Input function parameters
  172. ;    x = findgen(nx) # replicate(1.0, ny)    ;Create X and Y arrays
  173. ;    y = replicate(1.0, nx) # findgen(ny)
  174. ;    u = ((x-a(4))/a(2))^2 + ((y-a(5))/a(3))^2  ;Create ellipse
  175. ;    z = a(0) + a(1) * exp(-u/2)        ;to gaussian
  176. ;    z = z + randomn(seed, nx, ny)        ;Add random noise, SD = 1
  177. ;    yfit = gauss2dfit(z,b)            ;Fit the function, no rotation
  178. ;    print,'Should be:',string(a,format='(6f10.4)')  ;Report results..
  179. ;    print,'Is:      :',string(b(0:5),format='(6f10.4)')
  180. ;
  181. ; MODIFICATION HISTORY:
  182. ;    DMS, RSI, June, 1995.
  183. ;-
  184. ;
  185. on_error,2                      ;Return to caller if an error occurs
  186. s = size(z)
  187. if s[0] ne 2 then $
  188.     message, 'Z must have two dimensions'
  189. n = n_elements(z)
  190. nx = s[1]
  191. ny = s[2]
  192. np = n_params()
  193. if np lt 3 then x = findgen(nx)
  194. if np lt 4 then y = findgen(ny)
  195.  
  196. if nx ne n_elements(x) then $
  197.     message,'X array must have size equal to number of columns of Z'
  198. if ny ne n_elements(y) then $
  199.     message,'Y array must have size equal to number of rows of Z'
  200.  
  201. if keyword_set(neg) then q = MIN(SMOOTH(z,3), i) $
  202.     ELSE q = MAX(SMOOTH(z,3), i)    ;Dirty peak / valley finder
  203. ix = i mod nx
  204. iy = i / nx
  205. x0 = x[ix]
  206. y0 = y[iy]
  207.  
  208. xfit = gaussfit(x, z[*,iy], ax, NTERMS=4) ;Guess at params by taking slices
  209. yfit = gaussfit(y, z[ix,*], ay, NTERMS=4)
  210.  
  211. ; First guess, without XY term...
  212. a = [    (ax[3] + ay[3])/2., $        ;Constant
  213.     sqrt(abs(ax[0] * ay[0])), $    ;Exponential factor
  214.     ax[2], ay[2], ax[1], ay[1]]    ;Widths and centers
  215.  
  216. ;  If there's a tilt, add the XY term = 0
  217. if Keyword_set(tilt) then a = [a, 0.0]
  218.  
  219. ;************* print,'1st guess:',string(a,format='(8f10.4)')
  220. result = curvefit([nx, ny, x, y], reform(z, n, /OVERWRITE), $
  221.         replicate(1.,n), a, itmax=50, $
  222.         function_name = "GAUSS2_FUNCT", /NODERIVATIVE)
  223.  
  224. ; If we didn't already have an XY term, add it = 0.0
  225. if Keyword_set(tilt) eq 0 then a = [a, 0.0] $
  226. else a[6] = a[6] mod !pi        ;Reduce angle argument
  227. z= REFORM(z, nx, ny, /OVERWRITE)    ;Restore dimensions
  228. return, REFORM(result, nx, ny, /OVERWRITE)
  229. end
  230.